library(ggpubr)
library(anytime)
library(googlesheets4)
library(ggpmisc)
library(plyr)
library(dplyr)
library(lubridate)
library(tidyverse)
library(lubridate)
library(ggplot2)
library(scales)
library(zoo)
library(xts)
library(forecast)
library(googledrive)
library(streamMetabolizer)
library(readr)
library(here)
library(rstan)
library(imputeTS)
library(itsmr)
library(purrr)
library(naniar)
library(data.table)
library(grid)

#Metab 2019.1

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2019-05-20 04:13:57" & solar.time < "2019-07-01 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.266425 hours
save(mm.test.strt, file = here("Outputs", "strt.2019.1.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2019.1.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2019.1.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 2 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 2 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 2 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 2 rows containing non-finite values (stat_poly_line).
## Warning: Removed 2 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 2 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 2 rows containing non-finite values (stat_poly_line).
## Warning: Removed 2 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 2 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 2 rows containing non-finite values (stat_poly_line).
## Warning: Removed 2 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 2 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##     -0.6768       1.0619
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).

 rm(list = ls())

#Metab 2019.2

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2019-06-30 04:13:57" & solar.time < "2019-12-31 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.731593 hours
save(mm.test.strt, file = here("Outputs", "strt.2019.2.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2019.2.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2019.2.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 1 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##     -0.3249       1.0303
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

#Metab 2020.1

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2020-06-15 04:13:57" & solar.time < "2020-08-04 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.935815 hours
save(mm.test.strt, file = here("Outputs", "strt.2020.1.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2020.1.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2020.1.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 1 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##      -1.190        1.103
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

 rm(list = ls())

#Metab 2020.2

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2020-08-02 04:13:57" & solar.time < "2020-09-04 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.356035 hours
save(mm.test.strt, file = here("Outputs", "strt.2020.2.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2020.2.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2020.2.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##      -2.705        1.230
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))

 rm(list = ls())

#Metab 2020.3

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2020-09-02 04:13:57" & solar.time < "2020-12-04 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.394233 hours
save(mm.test.strt, file = here("Outputs", "strt.2020.3.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2020.3.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2020.3.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 1 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 1 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 1 rows containing non-finite values (stat_poly_line).
## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 1 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##     -0.3065       1.0220
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).
## Removed 1 rows containing missing values (geom_point).

 rm(list = ls())

#Metab 2021.1

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2021-05-05 04:13:57" & solar.time < "2021-06-20 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 1.81814 hours
save(mm.test.strt, file = here("Outputs", "strt.2021.1.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2021.1.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2021.1.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 3 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##      -1.092        1.088
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

 rm(list = ls())

#Metab 2021.1

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2021-05-05 04:13:57" & solar.time < "2021-06-23 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 2.200225 hours
save(mm.test.strt, file = here("Outputs", "strt.2021.1.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2021.1.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2021.1.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 3 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##      -1.072        1.087
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

 rm(list = ls())

#Metab 2021.2

bayes_name <- mm_name(type='bayes', pool_K600='binned', err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_name_AppEtAL <- "b_Kb_oipi_tr_plrckm.stan"
bayes_specs <- specs(bayes_name, 
                     burnin_steps=13000, saved_steps=2000, n_cores=4, n_chains = 4, GPP_daily_lower = 0, ER_daily_upper = 0
)
data.strt.mm.all <- read.csv(here("outputs", "stuart.comb.readyforMetab.new.clean.csv"))

data.strt.mm.all <- data.strt.mm.all %>% filter(solar.time >= "2021-06-22 04:13:57" & solar.time < "2021-12-01 04:13:57")

data.strt.mm.all <- data.strt.mm.all %>%  rename(DO.sat = DO.sat.EXO)
data.strt.mm.all$solar.time <- as.POSIXct(data.strt.mm.all$solar.time, tz = "UTC")
data.strt.mm.all <- data.strt.mm.all %>% select(solar.time, DO.obs, DO.sat, light, discharge, depth, temp.water)

startTime <- as.POSIXct(Sys.time())
mm.test.strt <- metab(bayes_specs, data=data.strt.mm.all)
endTime <- as.POSIXct(Sys.time())
# 
difftime(endTime, startTime)
## Time difference of 24.17826 mins
save(mm.test.strt, file = here("Outputs", "strt.2021.2.new.RData"))
# 
fit.strt <- get_fit(mm.test.strt)
fit.daily <- fit.strt$daily
write.csv(fit.daily, here("outputs", "strt.2021.2.new.full.csv"))

#Plot Input Data

model_data <- get_data(mm.test.strt)
write.csv(model_data, here("outputs", "strt.2021.2.new.model_data.csv"))

plot_DO_preds(mm.test.strt)

light.fig <- ggplot(data=model_data, aes(x=solar.time, y=light)) + geom_point(color = "chartreuse4") +theme(text = element_text(size=20))
q.fig<- ggplot(data=model_data, aes(x=solar.time, y=discharge)) + geom_point(color = "darksalmon")+theme(text = element_text(size=20))
temp.fig<- ggplot(data=model_data, aes(x=solar.time, y=temp.water)) + geom_point(color = "turquoise4")+theme(text = element_text(size=20))
depth.fig<- ggplot(data=model_data, aes(x=solar.time, y=depth)) + geom_point(color = "maroon")+theme(text = element_text(size=20))

grid.newpage()
grid.draw(rbind(ggplotGrob(light.fig), ggplotGrob(q.fig),ggplotGrob(depth.fig), ggplotGrob(temp.fig)))

gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4")
gpp
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
gpp <- ggplot(data=fit.strt$daily, aes(x=date, y=GPP_mean)) + geom_point(color = "chartreuse4") + geom_errorbar(aes(ymin=GPP_mean-GPP_sd, ymax=GPP_mean+GPP_sd), width=.2, position=position_dodge(0.05), color = "chartreuse4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
er <- ggplot(data = fit.daily, aes(x = date)) +geom_point(aes(y = ER_mean), color = "firebrick3") + geom_errorbar(aes(ymin=ER_mean-ER_sd, ymax=ER_mean+ER_sd), width=.2, position=position_dodge(0.05), color = "firebrick3")
er
## Warning: Removed 3 rows containing missing values (geom_point).

k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange")
k600
## Warning: Removed 3 rows containing missing values (geom_point).

#save with no axis for combined plot
k600 <- ggplot(data=fit.strt$daily, aes(x=date, y=K600_daily_mean)) + geom_point(color = "orange") + geom_errorbar(aes(ymin=K600_daily_mean-K600_daily_sd, ymax=K600_daily_mean+K600_daily_sd), width=.2, position=position_dodge(0.05), color = "orange") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )
rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
rhat
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

rhat <- ggplot(data=fit.strt$daily, aes(x=date)) + geom_point(aes(y=GPP_Rhat, colour = "GPP")) + geom_point(aes(y=as.numeric(ER_Rhat), colour = "ER")) + geom_point(aes(y=K600_daily_Rhat, colour = "K600")) +
  labs(y = "RHAT") + theme(legend.position="top") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() )+ geom_hline(yintercept=1.01, color = "dark blue")+ geom_hline(yintercept=1.1, color = "dark red", linetype = "dashed")
ggplot(data=fit.daily, aes(x=GPP_mean, y=ER_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=GPP_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq() 
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(data=fit.daily, aes(x=ER_mean, y=K600_daily_mean)) + geom_point()+ stat_poly_line() + stat_poly_eq()
## Warning: Removed 3 rows containing non-finite values (stat_poly_line).
## Warning: Removed 3 rows containing non-finite values (stat_poly_eq).
## Warning: Removed 3 rows containing missing values (geom_point).

mcmc <- get_mcmc(mm.test.strt)
k600mcmc <- rstan::traceplot(mcmc, pars='K600_daily', nrow=7)
k600mcmc

gppmcmc <- rstan::traceplot(mcmc, pars='GPP_daily', nrow=7)
gppmcmc

ermcmc <- rstan::traceplot(mcmc, pars='ER_daily', nrow=7)
ermcmc

#DO r2

strt.DO.pred <- predict_DO(mm.test.strt)
strt.DO.pred <- na.omit(strt.DO.pred)
lm(DO.obs~ DO.mod, data = strt.DO.pred)
## 
## Call:
## lm(formula = DO.obs ~ DO.mod, data = strt.DO.pred)
## 
## Coefficients:
## (Intercept)       DO.mod  
##      -3.680        1.336
strt.DO.pred$SM.date <- as.Date(strt.DO.pred$solar.time - 4*60*60)
strt.DO.pred$SM.date <- as.factor(strt.DO.pred$SM.date)
fit_model <- function(x, y) summary(lm(x ~ y))$adj.r.squared
test.run <- strt.DO.pred %>% group_by(SM.date) %>% summarise(adj.R2 = fit_model(DO.obs,DO.mod))
strt.DO.pred <- full_join(strt.DO.pred, test.run)
## Joining, by = "SM.date"
rsq <-  summary(lm(strt.DO.pred$DO.mod~strt.DO.pred$DO.obs))$r.squared

p1 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=DO.obs, colour = "Observed DO"), color = "darkcyan") + geom_line(aes(y = DO.mod, colour= "Modeled DO"), color = "blue4") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank() ) + labs(colour="Green", y = "DO (mg/L)") +  theme(legend.position="none") + ggtitle("Observed (points) and Modeled (lines)")
p2 <- ggplot(strt.DO.pred, aes(x = solar.time)) + geom_point(aes(y=adj.R2))
q.fig <- ggplot(data = data.strt.mm.all, aes(x = solar.time)) +geom_point(aes(y=discharge))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2)))

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(rhat),ggplotGrob(q.fig), ggplotGrob(k600),ggplotGrob(gpp), ggplotGrob(er)))
## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).
## Removed 3 rows containing missing values (geom_point).

 rm(list = ls())